library("ggplot2")

make_tstat_plot <- function(..., ylab = "") {
  ggplot(...) +
    geom_histogram(aes(y = after_stat(density)),
                   breaks = seq(0, 2 * qnorm(0.975), qnorm(0.975) / 18),
                   colour = "black", 
                   fill = "white",
                   linewidth = 0.65) +
    ylim(0, 2) +
    geom_vline(aes(xintercept = qnorm(0.975)),
               color = "red", 
               linewidth = 1,
               alpha = 0.6) +
    scale_x_continuous(breaks = seq(qnorm(0.975) / 2, 2 * qnorm(0.975), qnorm(0.975) / 2),
                       limits = c(0, 2 * qnorm(0.975)),
                       labels = c(0.98, 1.96, 2.94, 3.92)) +
    theme_classic() +
    theme(plot.title = element_text(hjust = 0.5, size = 10)) +
    labs(title = "",  x = "t-statistic", y = ylab)
}

save_tstat_plot <- function(path, plot) {
  ggsave(path, plot = plot, height = 4.1, width = 4.1, units = "in")
}

rd_papers <- read.csv("data/cl-rd-papers-meta.csv")
reanalysis_res <- read.csv("data/reanalysis-results.csv")

col_to_merge <- c("estimate_code", "bw_select_strat", "tstat", "weight")
reanalysis_res <- merge(
  reanalysis_res,
  subset(rd_papers, data_avail, col_to_merge),
  by = "estimate_code"
)
rm(col_to_merge)


# Figure 1: Distribution of t-statistics among published RD studies in political science

# Panel A
fig1a <- make_tstat_plot(
  subset(rd_papers, main_estimate),
  aes(x = tstat),
  ylab = "Density"
)
save_tstat_plot("output/figures/figure1a-tstat-orginal-just-one.pdf", fig1a)

# Panel B
fig1b <- make_tstat_plot(
  subset(rd_papers, disag_estimate),
  aes(x = tstat, weight = weight)
)
save_tstat_plot("output/figures/figure1b-tstat-orginal-all-weighted.pdf", fig1b)


# Figure 2: Distributions of t-statistics disaggregated by bandwidth selection procedure

# Panel A
fig2a <- make_tstat_plot(
  subset(rd_papers, disag_estimate & bw_select_strat != "arbitrary"),
  aes(x = tstat, weight = weight),
  ylab = "Density"
)
save_tstat_plot("output/figures/figure2a-tstat-automated.pdf", fig2a)

# Panel B
fig2b <- make_tstat_plot(
  subset(rd_papers, disag_estimate & bw_select_strat == "arbitrary"),
  aes(x = tstat, weight = weight)
)
save_tstat_plot("output/figures/figure2b-tstat-arbitrary.pdf", fig2b)


# Figure 3: Distributions of t-statistics among replicated studies by method of analysis

# Panel A
fig3a <- make_tstat_plot(
  subset(reanalysis_res, bw_select_strat != "cct mse optimal"),
  aes(x = tstat, weight = weight),
  ylab = "Density"
)
save_tstat_plot("output/figures/figure3a-tstat-reanalysis-original.pdf", fig3a)

# Panel B
fig3b <- make_tstat_plot(
  subset(reanalysis_res, bw_select_strat != "cct mse optimal"),
  aes(x = cct_conven_tstat, weight = weight)
)
save_tstat_plot("output/figures/figure3b-tstat-reanalysis-conven.pdf", fig3b)

# Panel C
fig3c <- make_tstat_plot(
  subset(reanalysis_res, bw_select_strat != "cct mse optimal"),
  aes(x = cct_robust_tstat, weight = weight)
)
save_tstat_plot("output/figures/figure3c-tstat-reanalysis-robust.pdf", fig3c)


# Figure 5: Distributions of t-statistics disaggregated by data source

# Panel A
fig5a <- make_tstat_plot(
  subset(rd_papers, disag_estimate & election_data),
  aes(x = tstat, weight = weight),
  ylab = "Density"
)
save_tstat_plot("output/figures/figure5a-tstat-election.pdf", fig5a)

# Panel B
fig5b <- make_tstat_plot(
  subset(rd_papers, disag_estimate & !election_data),
  aes(x = tstat, weight = weight)
)
save_tstat_plot("output/figures/figure5b-tstat-nonelection.pdf", fig5b)


# Figure S3: Distributions of t-statistics among all studies by method of analysis

# Panel A
figS3a <- make_tstat_plot(
  reanalysis_res,
  aes(x = tstat, weight = weight),
  ylab = "Density"
)
save_tstat_plot("output/figures/figureS3a-tstat-reanalysis-all-original.pdf", figS3a)

# Panel B
figS3b <- make_tstat_plot(
  reanalysis_res,
  aes(x = cct_conven_tstat, weight = weight)
)
save_tstat_plot("output/figures/figureS3b-tstat-reanalysis-all-conven.pdf", figS3b)

# Panel C
figS3c <- make_tstat_plot(
  reanalysis_res,
  aes(x = cct_robust_tstat, weight = weight)
)
save_tstat_plot("output/figures/figureS3c-tstat-reanalysis-all-robust.pdf", figS3c)


# Figure S5 & S6: RDHonest results

figS5 <- make_tstat_plot(
  subset(reanalysis_res, !is.na(rdhonest_tstat)),
  aes(x = tstat, weight = weight),
  ylab = "Density"
)
save_tstat_plot("output/figures/figureS5-tstat-rdhonest-original.pdf", figS5)

figS6 <- make_tstat_plot(
  subset(reanalysis_res, !is.na(rdhonest_tstat)),
  aes(x = rdhonest_tstat, weight = weight),
  ylab = "Density"
)
save_tstat_plot("output/figures/figureS6-tstat-rdhonest-rdhonest.pdf", figS6)
